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surface  geometry,  the  displacement  functions  and  elemental  degrees 
of  freedom,  and  the  modification  of  the  generalized  stiffness  method 
required  for  the  implementation  of  the  triangular  element  are 
described.  The  correlation  of  both  theoretical  and  experimental  results 
with  those  obtained  by  the  present  method  are  shown  along  with  ideal¬ 
izations  required  for  accurate  results.  Static  and  dynamic  solution 
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SECTION  I 
INTRODUCTION 

The  research  reported  in  this  paper  represents  a  portion  of  the  work  accomplished  under 
Contract  No.  F33615-67-C-1661,  sponsored  by  the  Flight  Dynamics  Laboratory,  Wright- 
Patterson  Air  Force  Base.  The  contract  work  has  as  its  goal  the  development,  application, 
and  evaluation  of  doubly-curved  finite  elements  in  the  dynamic  analysis  of  shells.  It  has  en¬ 
compassed  the  development  of  stiffness  and  mass  matrices  for  shell  elements  on  a  doubly 
curved  surface.  The  shell  geometry  is  limited  only  by  the  condition  that  the  surface  equations 
be  given  in  the  parametric  form 

x  =  f,  (a,  ,a2) 
y  =  f2(a,  ,a2) 

z  =f3(a.>  a2> 

in  which  f^,  f f^  are  known  functions  of  the  parameters  and  a^>  and  the  latter  con¬ 
stitute  an  orthogonal,  principal  curvilinear  coordinate  system  on  the  surface  of  the  shell. 

Derivation  of  the  shell  element  stiffness  matrices  is  based  on  the  strain  displacement 
equations  of  Novozhilov,  Reference  1,  while  the  elemental  mass  matrices  are  derived  by 
classical  means  from  the  kinetic  energy  associated  with  the  admissible  displacement  states 
as  described  by  Archer,  Reference  2. 

Two  basic  element  planforms  are  considered;  one  is  a  quadrilateral  bounded  by  principal 
coordinate  lines  (hence,  a  “rectangle”  in  the  orthogonal,  curvilinear  coordinate  system), 
and  the  other  is  a  triangle  with  arbitrarily  oriented  sides.  Supplementary  to  these  two  basic 
shell  elements  there  has  also  been  developed  an  apex  element,  which  is  a  degenerate  case 
of  the  quadrilateral.  It  is  obtained  by  modifying  the  quadrilateral  element  to  make  it  suitable 
for  use  at  the  apex  of  a  shell  of  revolution,  where  there  is  a  pole  in  the  surface  coordinate 
system.  A  stiffener  element  lying  on  a  meridonal  or  circumferential  coordinate  line  of  a  sur¬ 
face  of  revolution  has  also  been  developed.  This  element  has  stretching,  bending,  and  torsional 
properties,  and  its  elastic  axis  may  be  offset  from  the  middle  surface  of  the  shell. 

To  apply  the  finite  elements  just  described  to  the  analysis  of  actual  shell  problems,  a  com¬ 
puter  code  has  been  developed.  The  code  generates,  from  basic  input  data  describing  the 
geometry,  material,  and  loading  of  the  shell,  the  elemental  stiffness,  mass,  and  loads  matrices 
for  the  above  elements,  and  assembles  them  into  gross  matrices  representing  the  linear,  un¬ 
damped  equations  of  motion  of  the  complete  shell.  These  equations  are  solved  for  static 
deflections  of  the  shell  and  for  natural  vibration  modes  and  frequencies. 
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Evaluation  of  the  shell  finite  element  analysis  is  being  accomplished  by  correlating  the 
results  obtained  from  the  computer  code  described  above  for  several  specific  shell  response 
problems,  both  static  and  dynamic,  with  exact  theoretical  solutions,  experimental  data,  or 
other  numerical  solutions. 

The  purpose  of  this  paper  is  to  describe  briefly  the  methods  developed  and  to  demonstrate 
their  validity  by  presenting  the  correlation  studies  which  have  been  completed  at  this  time. 
The  description  of  the  methods  is  limited  to  what  are  considered  the  most  important  concepts, 
which  include  the  treatment  of  the  shell  surface  geometry,  the  displacement  functions  and 
elemental  degrees  of  freedom,  and  the  modifications  to  the  existing  finite  element  technology 
which  are  required  for  the  implementation  of  the  triangular  element.  Numerical  results  for 
shells  analyzed  using  stiffener  elements  and  triangular  shell  elements  are  not  available  at  the 
time  of  writing  this  paper;  therefore,  the  examples  presented  are  limited  to  unstiffened  shells 
and  to  the  use  of  quadrilateral  elements. 

More  detailed  descriptions  of  the  methods  developed,  including  detailed  element  deriva¬ 
tions,  will  appear  in  the  contract  final  report,  which  will  be  completed  later  this  year. 


187 


A  FFDL-TR- 68-150 


SECTION  II 

SHELL  ELEMENT  GEOMETRY 

The  computer  code  must  be  applicable  to  a  very  wide  range  of  doubly  or  singly  curved 
shells.  In  the  computation  process,  an  approximate  shell  geometry  must  be  defined  from  given 
input  data.  Because  of  the  great  variety  of  shell  shapes  and  the  rather  complex  forms  of  their 
geometrical  descriptions,  in  terms  of  a  set  of  interrelated  mathematical  functions,  this  geo¬ 
metrical  definition  is  a  difficult  task.  There  are  two  alternative  procedures  which  could  be 
used.  The  first  is  the  construction  of  an  assemblage  of  shell  elements  which  fits  together 
in  a  continuous  fashion  and  is  a  close  approximation  to  the  shape  of  the  actual  shell.  This 
almost  certainly  would  require  a  matching  of  nodal  locations  between  the  element  assemblage 
and  the  actual  shell.  It  should  also  include  a  matching  of  the  two  slope  components  at  each 
node.  Hopefully,  each  element  of  the  assemblage  could  have  a  relatively  simple  doubly  curved 
shape.  This  approach  is  extremely  difficult,  in  general,  for  the  reason  that  simple  element 
shapes  do  not  permit  sufficiently  continuous  assemblages  for  actual  shells  of  a  wide  range  of 
shapes.  The  second  procedure  is  to  approximate,  in  a  physical  way,  the  shell  itself.  The 
geometrical  functions  which  are  approximated  are  the  Lame  parameters,  A^,  A2>  and  A^. 
which  determine  length  and  angle  relationships  in  the  shell  middle  surface,  and  the  radii  of 
curvature,  R^,  R^,  and  the  twist,  R^,  which  determine  the  local  departure  of  the  shell 
from  planeness.  By  representing  these  functions  by  approximate  forms,  such  as  polynomials, 
in  the  shell  middle  surface  coordinate  system,  a  close  approximation  to  the  actual  shell  is 
obtained  in  a  numerical  sense.  However,  the  approximate  functions  do  not  represent  any  actual 
shell  at  all;  that  is,  they  do  not  necessarily  satisfy  the  surface  compatibility  equations  of 
Gauss  and  Codazzi. 

In  addition  to  the  six  quantities  listed  above,  their  derivatives  with  respect  to  the 
surface  coordinates,  as  well  as  the  shell  thickness,  t,  appear  as  products  in  various  com¬ 
binations  within  the  terms  of  the  strain  energy  integral.  To  evaluate  the  energy  accurately, 
each  such  term  should  be  approximated  individually,  starting  with  the  true  shell  functions, 
rather  than  by  taking  derivatives  and  forming  products  of  approximate  functions  already  de¬ 
fined. 

The  purpose  of  the  geometrical  representation  is  to  make  possible  the  evaluation  of  the 
strain  energy  integral  for  the  shell  elements.  It  is  believed  that  an  accurate  evaluation  of  the 
strain  energy  function  is  more  important  than  adherence  to  a  true  shell  surface,  composed 
of  elements  each  of  which  may  depart  significantly  in  local  geometry  from  the  actual  shell. 
For  this  reason,  the  second  approach  has  been  used. 
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All  derivations  and  calculations  have  been  based  on  the  use  of  a  principal  coordinate 
system  in  the  shell  middle  surface,  defined  by  the  variables  a,  and  a  In  this  case  A10 
an<3  Rj2  vanis^>  and  R^»  R2  are  the  principal  radii  of  curvature  of  the  mid-surface.  The  nec¬ 
essary  actual  shell  description  for  problem  solutions  consists  of  the  functions  (either  in 
mathematical  or  dense  tabular  form)  Ar  Ag,  R  ,  R2,  and  t  in  terms  of  the  principal  coordinates 
of  the  surface.  In  addition,  it  is  necessary  to  have  a  description  of  the  principal  coordinate 
system  in  terms  of  a  base  rectangular  cartesian  reference  frame,  in  order  to  determine  the 
locations  of  loadings,  nodes,  etc.,  in  terms  of  and  a  The  approximations  to  the  terms 
of  the  strain  energy  integral,  in  terms  of  Qt  ^  and  a  are  second  degree  polynomials 
determined  by  matching  the  individual  function  values  at  points  on  the  boundaries  and  in  the 
interior  of  the  elements. 


SECTION  m 

SHELL  ELEMENT  DISPLACEMENT  FUNCTIONS 

The  element  displacement  functions  were  chosen  on  the  basis  of  an  accumulation  of  re¬ 
quirements  which  has  gradually  taken  shape  as  a  result  of  the  work  of  a  number  of  re¬ 
searchers,  It  appears  worthwhile  to  discuss  this  accumulation  of  knowledge  here,  since  such 
a  summary  has  not  been  given  previously, 

1.  The  rigid  body  motion  of  an  element  must  be  represented,  at  least  approximately, 
in  the  displacement  functions.  This  item  has  been  the  subject  of  some  controversy.  Jones 
and  Strome,  (Reference  3),  included  rigid  motions  explicitly  in  their  doubly  curved  shell 
element  and  found  a  gain  in  accuracy  to  occur.  Cantin  and  Clough,  (Reference  4)  found  similar 
results  in  their  cylindrical  shell  element.  On  the  other  hand,  Bogner,  Fox,  and  Schmidt, 
(Reference  5)  obtained  excellent  results  with  a  cylindrical  shell  element  without  explicit 
inclusion  of  rigid  motions.  Haisler  and  Stricklin,  (Reference  6)  also  presented  data  supporting 
the  omission  of  rigid  modes.  The  conclusion  which  appears  to  derive  from  all  work  to  date  is 
that  rigid  motions  must  be  represented,  but  that  this  can  be  done  adequately  by  means  of 
polynomial  displacement  functions  of  sufficiently  high  degree. 
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2.  The  functions  u,  v,  and  w  should  all  be  represented  with  equal  accuracy,  i,e.,  by 
polynomials  of  equal  degree.  The  basis  of  this  is  also  contained  in  the  work  of  a  number  of 
authors.  Jones  and  Strome  found  that  the  explicit  inclusion  of  the  rigid  body  displacements  was 
insufficient  to  provide  either  a  high  level  of  accuracy  for  large  element  sizes  or  a  rapid 
convergence  with  decreasing  element  sizes.  In  order  to  obtain  such  results,  they  found  it 
necessary  to  add  two  other  displacement  forms.  These  amounted  to  explicit  inclusion  of  uniform 
hoop  expansion  and  rotation  of  an  axisymmetric  element  about  its  own  circumferential  central 
circle,  similar  in  nature  to  torsion  of  a  ring.  These  functions  had  the  effect  of  increasing  the 
degree  of  the  displacement  polynomials,  primarily  the  tangential  displacement  function,  which 
initially  had  been  of  the  first  degree  in  the  arc  length  variable.  With  these  functions,  good 
accuracy  has  been  obtained  for  elements  of  very  large  size.  The  Haisler  and  Stricklin  element 
has  the  same  polynomial  functions  as  the  initial  polynomial  of  Jones  and  Strome.  It  has  shown 
relatively  slow  convergence  and  the  need  for  small  elements  to  obtain  good  accuracy.  The 
Cantin  and  Clough  displacements  and  those  of  Bogner,  Fox,  and  Schmidt  differ  only  in  the  use, 
in  the  former,  of  a  bilinear  polynomial  for  u  and  v,  compared  to  a  bicubic  polynomial  in  the 
latter.  The  Bogner,  et.  al. ,  element  shows  much  faster  convergence  and  better  accuracy  for 
large  elements  than  that  of  Cantin  and  Clough,  Similar  conclusions  arise  from  the  results  of 
Greene,  et.  al. ,  (Reference  7).  Based  on  these  results,  the  need  for  equally  competent 
forms  for  all  of  u,  v,  and  w  appears  to  be  justified. 


3.  The  displacement  functions  for  u,  v,  and  w  should  be  at  least  of  the  competence  of  the 
bicubic  form 


4  4 

£  £  ^mn  ^  9n  ^ 

m  =  i  n=i 


where  fm  and  gn  are  descriptively  called  beam  functions.  That  is,  each  is  a  cubic  polynomial 
and  defines  the  displacement  within  a  line  element  in  terms  of  the  displacement  and  slope  at 
its  ends.  They  provide  the  ability  of  each  function  of  a1,  i.e.,  fm  (  a  jb  to  build  up  and  die 
out  in  the  a  ^  direction  in  one  element  span,  and  vice  versa. 


Several  investigators  have  independently  discovered  the  bicubic  functions,  but  most  have 
applied  it  incorrectly  or  in  an  incomplete  form.  The  correct  form  and  usage  is  given  for 
plate  elements  in  Reference  8  and  for  cylindrical  shell  elements  in  Reference  5.  Results  of 
good  accuracy  were  obtained  in  this  work.  Based  on  this  and  on  similar  results  obtained  by 
Greene,  et.  al. ,  the  bicubic  functions  have  been  adopted  for  the  present  work  on  doubly 
curved  shells. 
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4.  The  continuity  of  displacement  and  slope  must  be  enforced,  at  least  approximately,  at 
the  interelement  boundaries.  This  requirement  has  been  discussed  by  a  number  of  investigators. 
It  was  covered  most  completely  in  the  work  of  Greene,  et.  al.,  where  it  was  found  that  exact 
continuity  need  not  be  enforced  to  obtain  good  accuracy.  Instead,  with  competent  (reasonably 
high  degree)  displacement  functions,  good  accuracy  can  be  obtained  with  continuity  enforced 
only  on  the  lower  degree  polynomial  components  of  the  displacements  and  slopes  at  the  inter¬ 
element  boundaries.  The  bicubic  forms  provide  the  ability  to  obtain  exact  continuity  for  shell 
elements  whose  corners  are  right  angles,  provided  that  the  radius  of  curvature  is  continuous 
across  interelement  boundaries.  An  ingenious  set  of  functions  which  satisfies  explicitly  all 
continuity  requirements  for  a  flat  triangular  plate  element  is  given  in  References  9  and  10. 
Its  extension  to  a  doubly  curved  shell  element  appears  at  the  present  to  be  quite  complicated, 
and  therefore,  for  triangular  elements  it  was  decided  to  apply  the  generalized  finite  element 
approach,  using  Lagrange  multiplier  functions  to  enforce  approximate  continuity.  This  was 
described  by  Jones,  (Reference  11)  and  applied  by  Greene,  et.  al.,  (Reference  7),  An  ex¬ 
tension  of  this  approach  was  given  by  McLay,  (Reference  12)  which  greatly  increases  its 
applicability  to  the  practical  solution  of  problems. 

The  displacement  functions  used  for  each  of  u,  v,  and  w,  are  cubic  polynomials  in  two 
variables, 

3  3  _  _ 

I  I  BmnC  £ 

m=o  n  =  o 

where  £  and  £  are,  respectively,  dimensionless  variables  representing  a  and  a  0  and 
having  unit  change  over  the  total  dimensions  of  the  element.  The  element  matrix  derivations 
are  carried  out  in  the  terms  of  the  coefficients  B  ,  and  are  then  transformed  by  invariance 
to  the  more  physically  convenient  coefficients,  A  ,  of  the  previously  given  displacement 
form.  This  procedure  makes  it  convenient  to  carry  out  the  derivation  of  the  stiffness  matrix 
in  closed  form,  a  single  algebraic  expression  being  obtained  for  each  of  the  six  distinct 
partitions  of  the  matrix  when  it  is  partitioned  according  to  u,  v,  and  w  displacement  fields. 
The  transformed  degrees  of  freedom,  A  ,  include  at  each  node  of  a  quadrilateral  element 
for  each  of  u,  v,  and  w,  its  value,  its  derivative  with  respect  to  each  of  £  and  £  ,  and  its 
second  cross  derivative  with  respect  to  £  and  £  .  These  total  48  degrees  of  freedom  per 
element.  By  matching  between  elements  common  to  a  node  all  of  the  12  degrees  of  freedom 
per  node,  complete  continuity  of  displacements  and  their  first  derivatives  is  obtained  at  all 
interelement  boundaries.  The  continuity  of  the  derivatives  of  u  and  v  is  not  required,  how¬ 
ever,  for  potential  energy  formulations  such  as  the  finite  element  method.  For  this  reason, 
and  because  the  present  work  is  aimed  at  stiffened  shells,  in  which  the  derivatives  of  u  and 
v  should  actually  experience  jumps  at  element  boundaries  where  stiffeners  are  attached. 
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these  derivatives  are  not  made  continuous.  Omitting  this  requirement  is  accomplished  by 

-.2  %2 

omitting  the  nodal  matching  of  $  “  and  P  r  between  neighboring  elements.  These 

OQ  oC  dQ 

degrees  of  freedom  are  handled  by  eliminating  them  at  the  elemental  level.  They  become 
dependent  freedoms,  determined  by  the  values  of  the  remaining  40  freedoms  of  the  element. 
All  of  the  latter  freedoms  are  merged,  that  is,  matched  between  elements,  at  the  nodes.  The 
method  of  Guyan,  (Reference  13)  is  used  for  the  corresponding  reduction  of  the  mass  matrix. 

The  difference  between  the  displacement  functions  used  here  and  those  of  Reference  5 
lies  in  the  omission  here  of  the  requirement  of  continuity  of  the  u  and  v  derivatives  and  in 
the  transformation  to  the  more  physically  appropriate  freedoms  of  rotations  and  twist  of  the 
shell  atthe  nodes  rather  than  simply  the  derivatives  of  u,  v,  and  w.  The  latter  step  is  necessary 
to  preserve  the  continuity  of  slope  when  the  curvatures  of  the  shell  may  be  discontinuous  across 
interelement  boundaries.  It  is  also  helpful  in  merging  stiffener  element  freedoms  with  those  of 
the  shell. 

The  above  discussion  applies  primarily  to  rectangular  shell  elements,  which  are  expected 
to  form  the  bulk  of  the  elements  inproblems.  Triangular  elements  will  be  used  where  needed, 
such  as  around  cutouts  or  at  boundaries  inclined  to  the  coordinate  lines.  They  will  be  used 
only  in  conjunction  with  a  primary  system  of  rectangular  elements  and,  consequently,  it  is 
most  convenient  to  use  the  same  displacement  functions  for  the  triangles  as  for  the  rectangles. 
This  permits  exact  enforcement  of  continuity  requirements  where  a  triangle  joins  a  rectangle 
and  approximate  enforcement  of  continuity,  through  the  generalized  method,  where  two  triangles 
join. 


Though  numerical  results  for  the  triangular  elements  are  not  included  in  this  paper,  a 
brief  discussion  of  the  approximate  continuity  enforcement  appears  to  be  pertinent.  It  is  given 
in  the  next  section. 
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SECTION  IV 

GENERALIZED  FINITE  ELEMENT  TECHNIQUE 


The  triangular  elements  necessarily  have  two  corners  where  the  sides  do  not  meet  at 
right  angles.  Because  of  this,  it  is  very  difficult,  perhaps  impossible,  for  doubly  curved 
elements  to  maintain  the  required  continuities  of  both  displacement  and  slope  across  inter¬ 
element  boundaries.  In  addition,  the  triangular  elements,  having  one  less  node  than  the  quadri¬ 
laterals,  do  not  provide  through  nodal  quantities  a  complete  definition  of  the  bicubic  displace¬ 
ment  forms.  The  generalized  finite  element  technique  provides  a  means  for  enforcing  approx¬ 
imate  continuity  along  the  sides  of  adjacent  elements,  where  exact  continuity  is  impossible  or 
cannot  be  enforced  through  the  matching  of  nodal  quantities.  The  method  involves  formulating 
the  finite  element  equations  without  including  explicit  continuity  enforcement,  and  subsequently 
imposing  continuity  through  the  use  of  Lagrange  multiplier  functions  and  integral  constraint 
equations.  The  constraint  equations  are  solved  along  with  the  usual  finite  element  equations. 
The  integral  constraints  can  provide  exact  continuity  if  this  is  possible;  otherwise  approximate 
continuity  is  obtained. 

The  generalized  technique  has  the  disadvantage  that  it  populates  almost  completely  the 
mass  and  stiffness  matrices.  Since  the  sparseness  of  these  matrices  in  usual  finite  element 
formulations  is  a  great  advantage  in  computing  ease  and  economy,  the  generalized  technique 
is  not  readily  applicable  in  its  originally  given  form,  Reference  11.  A  modification  of  the 
method  has  been  developed  by  McLay,  (Reference  12)  which  eliminates  this  difficulty.  It 
creates  a  set  of  imaginary  lines  which  lie  between  the  sides  of  adjacent  elements  where 
Lagrange  multipliers  are  to  be  used  to  enforce  continuity.  These  lines,  called  buffer  lines, 
are  given  degrees  of  freedom  which  are  defined  by  rotations  and  displacements  of  their  end 
points.  The  deformations  and  displacements  of  the  buffer  lines  aredefined  by  these  endpoints, 
or  nodal,  degrees  of  freedom.  These  nodal  values  become  the  actual  freedoms  of  the  finite 
element  formulation.  The  element  displacement  coefficients  are  matched  to  the  buffer  line 
nodal  freedoms  through  the  use  of  Lagrange  multipliers  and  thereby  eliminated  from  the 
calculations.  The  resulting  formulation  has  loosely  coupled  mass  and  stiffness  matrices  as  in 
the  usual  finite  element  formulation.  It  provides,  approximately  or  exactly,  depending  on  the 
item  in  question,  the  required  continuities  of  the  finite  element  method. 
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SECTION  V 

DISCUSSION  OF  NUMERICAL  RESULTS 


In  order  to  demonstrate  the  capability  of  the  shell  elements  to  model  accurately  shell¬ 
like  structures,  correlations  with  analytical  and  experimental  data  are  given  for  several 
static  and  dynamic  problems. 

HEMISPHERICAL  SHELL 

The  first  problem  considered  involves  edge  bending  behavior  of  a  shell  of  revolution. 
A  hemispherical  shell  was  subjected  to  an  edge  radial  shear  load  as  shown  in  Figure  1. 
For  the  particular  geometry  chosen,  the  region  of  stress  and  deformation  is  confined  to  a 
15-degree  segment  near  the  boundary.  To  cover  adequately  this  region,  a  30-degree  segment 
was  used  in  the  analysis.  Figure  1  shows  the  normal  displacement  and  rotation  respectively 
for  three  different  finite  element  idealizations.  The  idealizations  are:  three  elements  at  10 
degrees  each,  six  elements  at  five  degrees  each,  and  12  elements  at  two  and  one  half  degrees 
each.  The  figure  shows  that  all  three  idealizations  yield  results  of  essentially  equivalent 
accuracy.  The  results  at  the  base  of  the  shell  are  compared  with  the  theoretical  solutions 
given  in  References  1  and  14,  Excellent  accuracy  is  obtained,  as  seen  in  Table  1. 


TABLE  1 

PERCENT  ERRORS  AT  BASE  OF  SHELL  -  SHELL  COMPUTER 
PROGRAM  RESULTS  COMPARED  TO  THEORETICAL  SOLUTION 


IDEALIZATIONS 

NORMAL  DISPLACEMENT 

SLOPE  CHANGE 

3  Elements  at  10° 

-4.95  % 

-0.16  % 

6  Elements  at  5° 

-1.98  % 

-0.30  % 

12  Elements  at  2-1/  2° 

-0. 22  % 

-0.12  % 

ELLIPSOIDAL  SHELL 

Previous  results  with  curved  finite  element  shell  solutions,  Jones  and  Strome,  {Ref¬ 
erence  3)  have  shown  that  accuracy  may  deteriorate  seriously  for  shell  shapes  other  than 
spherical,  particularly  for  large  elements  or  inadequate  displacement  functions.  In  order  to 
demonstrate  the  applicability  of  the  curved  shell  elements  for  such  shells,  a  2:1  ellipsoid  is 
considered.  Figure  2  shows  comparative  results  for  normal  displacement  and  rotation  from 
membrane  theory,  from  the  curved  shell  of  revolution  element  of  Reference  3,  and  from  two 
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idealizations  of  the  quadrilateral  shell  elements  discussed  in  this  paper.  The  ellipsoidal  shell 
is  loaded  by  internal  pressure  and  allowed  to  rotate  at  the  supports.  The  shell  exhibits 
predominantly  membrane  behavior.  However,  the  bending  theory  is  required  to  compute  the 
rotation  of  the  meridional  tangent.  Figure  2  shows  the  excellent  agreement  of  all  solutions 
for  normal  displacement,  and  the  inadequacy  of  the  membrane  theory  and  the  eight-element 
quadrilateral  shell  idealization  for  predicting  the  rotations.  The  shell  of  revolution  element 
of  Reference  3  and  the  15-element  quadrilateral  element  solutions  are  essentially  identical 
and  yield  excellent  correlation  at  the  boundary  with  the  available  analytical  data.  The  latter 
are  given  in  Reference  15.  Comparison  of  finite  element  and  exact  data  is  given  in  Table  2. 


TABLE  2 

NORMAL  DISPLACEMENT  AND  ROTATION  AT  BASE  OF 
INTERNALLY  PRESSURIZED  ELLIPSOIDAL  SHELL 


Theory 

Ref.  15 

Membrane 

Theory 

Shell  of 

Revolution 

Ref.  3 

Finite  Element  Solution 

8  Elements 

15  Elements 

w  ~ 

inches 

-.24048 

-.230 

24048 

-. 24675 

-.24088 

/3  s  — 

radians 

4.4024xl0~4 

0 

4.4341xl0-4 

7„4940xl0~4 

4.4217X10"4 

VIBRATION  OF  A  CYLINDER 

The  natural  frequencies  and  modes  of  a  simply  supported  circular  cylinder  without 
axial  restraint  have  been  studied  using  the  quadrilateral  curved  shell  element.  The  con¬ 
figuration  of  the  cylinder  analyzed  is  given  in  Figure  3.  Theoretical  solutions  for  the  natural 
frequencies  of  this  cylinder  presented  by  Voss  in  Reference  16  are  compared  with  the  re¬ 
sults  from  the  finite  element  analysis. 

In  the  first  study  a  single  mode  was  examined.  This  was  done  by  representing  only  one 
quarter-wave  segment  of  the  shell  in  the  finite  element  analysis  and  imposing  kinematic 
boundary  conditions  at  the  edges  of  this  segment  which  are  compatible  with  the  symmetry  or 
antisymmetry  of  the  structural  response  as  indicated  in  Figure  3. 
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t=  0.04  in.  p  =0.  260X|0“3|b-sec//’in4 

Figure  3.  Analyzing  a  Specific  vibration  Mode  by  Using  One 
Quarter- Wave  Segment  of  Shell 


The  purpose  of  this  study  is  to  examine  the  effects  of  different  element  size-to-wave- 
length  ratios  and  different  ways  of  distributing  the  mass  upon  the  accuracy  of  the  analysis. 
The  mode  examined  is  the  N  =  8,  M  =  1  mode  which  has  eight  full  waves  around  the  circum¬ 
ference  of  the  cylinder  and  one-half  wave  down  the  length  as  indicated  in  Figure  3.  This  mode 
was  selected  because  it  is  the  fundamental  (i.e.,  lowest  frequency)  mode  for  this  cylinder  thus 
assuring  that  it  would  be  obtained  first  inthe  eigenvalue  -  eigenvector  solution  of  the  equations 
of  motion  rather  than  some  higher  order  mode  having  the  same  boundary  conditions.  The 
three  different  idealizations  ofthe  quarter- wave  segment  of  the  shell  which  were  analyzed  con¬ 
sisted  of  one,  four,  and  nine  elements  as  shown  in  Figure  4,  which  also  shows  the  boundary 
conditions  imposed  at  the  edges  of  the  segment. 
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Figure  4.  Idealization  of  Quarter-Wave  Segment  of  Shell  Showing  Boundary  Conditions 
Arrangement  of  Elements,  and  Node  Point  Identification 
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Two  different  methods  of  distributing  the  mass  are  also  included  in  this  study.  The 
first  of  these  is  the  generalized  mass  derived  by  classical  means  from  the  kinetic  energy 
associated  with  the  admissible  displacement  states  (i.e, ,  those  assumed  in  the  derivation  of 
the  stiffness  matrix).  The  mass  matrix  thus  obtained  is  the  same  as  the  consistent  mass 
matrix  derived  by  Archer,  Reference  2,  and  has  the  same  characteristic  population  and 
coupling  as  the  stiffness  matrix.  The  second  method  of  mass  distribution  studied  is  the  lumped 
mass  distribution  in  which  the  total  structural  mass  is  divided  uniformly  into  point  masses 
placed  at  all  nodes  where  translational  freedoms  are  retained.  The  size  of  the  eigenvalue 
problem  may  be  reduced  by  eliminating  selected  degrees  of  freedom,  retaining  those  which 
are  considered  to  be  most  important  in  representing  the  modal  response.  This  reduction  is 
accomplished  by  the  method  described  by  Guyan,  (Reference  13).  The  effect  of  different  de¬ 
grees  of  reduction  of  the  original  structural  freedoms  was  also  examined. 

Results  of  this  study  are  presented  in  Table  3.  Frequencies  obtained  using  both  the 
generalized  mass  and  lumped  mass  distributions  are  to  be  compared  with  the  exact  solution 
of  527  radians  per  second.  The  number  of  elements  and  the  number  of  total  structural  degrees 
of  freedom  (before  reduction)  are  indicated,  as  well  as  the  number  and  identification  of  degrees 
of  freedom  retained  (after  reduction)  in  the  eigenvalue  solution. 

Examine  first  the  results  obtained  with  the  generalized  mass  distribution.  For  one  element, 
with  a  total  of  fourteen  structural  degrees  of  freedom,  very  little  change  is  seen  in  the  com¬ 
puted  frequency  when  the  number  of  retained  freedoms  is  reduced  from  fourteen  to  five  con¬ 
sisting  only  of  the  physical  translations  and  rotations  at  the  node  points,  both  results  being 
approximately  2  1/2%  high.  As  further  reductions  are  made,  down  to  three  degrees  of  freedom 
consisting  of  the  translations  at  the  nodes  only,  and  finally  down  to  one  degree  of  freedom 
consisting  of  the  normal  component  of  translation  only,  the  error  increases  to  approximately 
3  1/2  and  5%  respectively.  In  the  four-element  idealization,  the  original  48  structural  degrees 
of  freedom  are  reduced  to  12,  retaining  all  translational  freedoms,  and  to  four,  retaining  only 
the  normal  translational  freedoms,  yielding  an  identical  result  which  is  less  than  1/2%  in 
error.  The  nine-element  idealization, containing  102  structural  degrees  of  freedom  originally, 
is  reduced  similarly  to  27,  retaining  all  translational  freedoms,  and  then  to  nine,  retaining  only 
the  normal  translational  freedoms,  both  cases  giving  the  exact  solution  of  527  radians  per 
second.  Finally,  the  four-element  and  nine- element  idealizations  were  reduced  to  one  and 
four  degrees  of  freedom  respectively,  retaining  only  the  normal  translational  freedoms  at 
alternate  nodes  (the  subscripts  on  w  in  Table  3  refer  to  the  node  numbers  indicated  in 
Figure  4  which  were  retained).  A  sharp  increase  in  the  error  is  noted  for  both  idealizations 
after  this  last  reduction  step.  It  is  apparent  that  for  the  same  size  eigenvalue  problem  a  more 
accurate  and  less  costly  result  is  obtained  by  using  the  next  coarser  structural  idealization 
and  retaining  all  w  freedoms. 
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TABLE  3 


FREQUENCIES  CALCULATED  FOR  CYLINDRICAL  SHELL 
MODE  OF  FIGURE  3  FROM  VARIOUS  STRUCTURAL 
IDEALIZATIONS  AND  MASS  DISTRIBUTIONS 
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Examining  next  the  results  obtained  from  the  lumped  mass  distributions  given  in  Table  3, 
it  is  apparent  that  good  results  were  obtained  in  all  cases.  This  surprising  result  is  due  to  the 
fortuitous  circumstance  or  compensating  errors,  the  inherent  overstiffness  of  the  stiffness 
matrix  being  offset  by  the  apparent  increase  in  flexibility  caused  by  lumping  the  masses  and 
thereby  concentrating  the  inertial  forces. 

The  same  cylinder,  Figure  3,  has  also  been  analyzed  by  taking  a  full  quarter  of  the  shell, 
i.e.,  a  90-degree  sector,  in  an  attempt  to  obtain  all  the  lower  frequency  modes  in  their  correct 
order.  Thirty-two  elements  were  used  over  the  quarter  of  the  shell,  eight  circumferentially 
and  four  axially.  Antisymmetric  boundary  conditions  were  prescribed  on  the  two  edges  of  the 
90-degree  sector.  This  eliminates  the  axial  rigid  body  mode  of  the  cylinder  and  allows  us  to 
obtain  the  modes  having  even  numbers  of  waves  circumferentially.  The  results  obtained  are 
plotted  in  Figure  5  as  frequency  versus  the  number  of  circumferential  waves,  N,  for  different 
values  of  M,  the  number  of  axial  half-waves,  as  indicated  by  the  three  distinct  dashed  curves. 
The  corresponding  theoretical  solutions  from  Reference  16,  are  shown  as  solid  curves.  This 
figure  shows  graphically  the  decrease  in  accuracy  as  the  element- size-to- wavelength  ratio 
increases  in  both  the  axial  and  circumferential  directions  as  indicated  by  increasing  values  of 
M  and  N  respectively.  Numerically,  this  error  increases  from  less  than  one-half  percent  at 
M  =  1,  N  =  4,  to  approximately  18  percent  at  M  =  3,  N  =  16.  The  16  modes  obtained  from  the 
finite  element  analysis  generally  appear  in  the  correct  sequence  of  ascending  frequency, 
although  some  are  interchanged  where  a  higher  order  mode  has  a  slightly  lower  frequency  than 
a  lower  order  mode.  In  these  cases,  the  increased  stiffness  of  the  finite  element  model  due  to 
the  larger  element- size-to- wavelength  ratio  has  caused  a  greater  increase  in  the  error  of  the 
finite  element  analysis  between  the  two  modes  than  the  theoretical  difference  in  frequency. 
For  example,  in  the  case  of  M  =  l,  the  N  =  6  and  the  N  =  10  modes  are  interchanged  as  are  the 
N  =  4  and  the  N  =  14  modes. 


The  freedoms  retained  in  the  quarter  shell  analysis  are  the  normal  displacement  com¬ 
ponent,  w,  and  the  rotation,  /3 q,  at  all  nodes.  Thus,  approximately  20%  of  the  original  structural 
degrees  of  freedom  were  retained  as  mass  points.  The  generalized  mass  distribution  described 
previously  was  used.  The  reason  for  retaining  in  this  study  was  to  better  represent  the 
higher  modes  having  N  greater  than  8,  thus  having  circumferential  element-size-to-wavelength 
ratios  greater  than  1:4,  as  in  the  coarsest  idealization  used  in  the  quarter-wave  study. 


To  compare  the  finite  element  vibration  analysis  with  experimental  results,  tests  of  an 
unstiffened  cylinder  with  fixed-free  boundary  conditions,  as  reported  in  Reference  17  were 
selected.  The  specific  model  is  that  referred  to  as  Model  1  in  the  report  and  consisted  of  a 
circular  cylinder  48  inches  long,  20  inches  in  diameter,  and  0.03  inches  thick,  fully  clamped 
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Figure  5.  Natural  Frequencies  of  a  Circular  Cylinder  Calculated  by 
Finite  Element  Analyses  Compared  With  Theory 
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Figure  6.  Correlation  of  Experimental  Versus  Analytical  Results  for 
Vibration  of  a  Fixed-Free  Cylindrical  Shell 


at  one  end  and  free  at  the  other.  The  finite  element  idealization  consisted  of  24  elements 
over  one  quarter  of  the  shell,  four  axially  and  six  circumferentially,  with  symmetric  boundary 
conditions  at  the  edges  of  the  quadrant  of  the  cylinder.  The  generalized  mass  distribution 
with  w  and  $q  freedoms  retained  at  each  node  was  used.  The  experimental  results  from 
Reference  17  and  the  first  six  modes  obtained  in  the  finite  element  analysis  are  shown 
in  Figure  6.  Note  that  here  the  definition  of  M  and  N  differs  from  that  given  previously,  M 
being  the  number  of  nodes  on  an  axial  line  (generator)  of  the  cylinder  and  N  being  the  number  of 
nodes  on  a  circumferential  line.  Nodes  are  understood  here  to  be  points  of  zero  normal 
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displacement,  w.  With  the  boundary  conditions  used  only  modes  having  even  numbers  ofwaves 
circumferentially  are  obtained,  thus  analytical  results  are  plotted  only  at  values  of  N  which 
are  multiples  of  four. 

The  correlation  is  quite  good;  the  increase  in  error  with  increase  of  M  and  N  noted  in 
the  correlation  with  theory  is  not  apparent  here.  This  may  be  due  in  part  to  the  fully  clamped 
boundary  conditions  assumed  in  the  analysis  at  the  fixed  end,  which  probably  differ  from  those 
actually  achieved  in  the  experiment.  This  might  very  well  influence  the  experimental  frequencies 
of  the  lower  order  modes. 

SPHERICAL  CAP 

Axi symmetrical  vibration  modes  were  computed  for  a<£=  60  degrees  spherical  cap  with 
hinged-fixed  boundaries  (i.e.,  zero  displacement  and  moment).  The  shell  and  its  properties 
are  shown  in  Figure  7,  Calculations  were  performed  for  two  idealizations  and  frequencies 
and  mode  shapes  are  compared  with  theory  in  Figures  8  to  11. 

The  spherical  cap  was  idealized  by  six  and  ten  elements  located  along  the  meridian.  Since 
the  finite  element  solutions  are  compared  only  with  theoretical  axisymmetric  modes,  the  shell 
was  idealized  by  a  15-degree  slice  with  proper  boundary  conditions  along  the  meridional  lines 
to  enforce  axisymmetric  modal  patterns.  The  six  and  ten  element  idealization  contained 
47  and  97  freedoms,  respectively,  with  the  stiffness  and  mass  matrices  reduced  to  11  and  19 
freedoms,  respectively,  for  the  eigensolution.  Only  the  normal  displacement  w  was  retained 
as  a  freedom  in  the  eigensolution  with  the  reduction  of  mass  matrix  calculation  based  on 
Guyan’s  method,  (Reference  13). 

Figures  8  through  11  give  the  first  four  theoretical  mode  shapes  and  frequencies  and  the 
finite  element  solutions  for  the  two  idealizations.  The  first  and  second  mode  shapes  are  as 
expected  with  one  and  two  crossings  of  the  meridional  line  respectively.  The  third  mode  is 
peculiar  in  that  it  does  not  cross  the  meridional  line  (i.e.,  no  interior  nodes).  This  mode  is 
missing  in  the  theoretical  solution  given  by  Kalnins,  (Reference  18),  but  was  later  determined 
as  a  proper  mode  shape  by  Cohen,  (Reference  19).  The  fourth  axisymmetric  mode  shape  con¬ 
tains  three  interior  nodes  and  corresponds  to  Kalnins’ s  third  and  Cohen’s  fourth  mode 
respectively.  Although  the  ten-element  idealization  gave  excellent  results  for  all  modes,  the 
six-element  idealization  experienced  difficulty  in  determining  the  third  mode  shape  accurately 
and  was  in  considerable  error  on  the  fourth  mode.  The  excessive  stiffening  that  occurs  in  the 
six-element  idealization  for  the  higher  modes  reflects  the  inadequacy  of  this  idealization  to 
represent  the  waveforms. 
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Figure  7.  Geometry,  Material  Properties,  and  Idealizations  -  Spherical  Cap 


Figure  8.  First  Axisymmetric  Mode  Shape  -  60°  Spherical  Cap 
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Figure  9.  Second  Axisymmetric  Mode  Shape  -  60°  Spherical  Cap 


a)  -rad/sec 


Figure  10.  Third  Axisymmetric  Mode  Shape  -  60°  Spherical  Cap 
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Figure  11.  Fourth  Axisymmetric  Mode  Shape  -  60°  Spherical  Cap 
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SECTION  V 
CONCLUSIONS 

A  capability  for  the  dynamic  analysis  of  doubly  curved  shells  by  the  finite  element  method 
has  been  developed.  The  derivation  of  the  shell  elements  and  their  implementation  have  been 
discussed  in  some  detail,  with  particular  emphasis  on  the  choice  of  displacement  functions, 
the  means  of  approximating  the  shell  geometry,  and  the  extension  of  the  generalized  finite 
element  method  which  is  being  used  in  applying  the  methods  to  a  triangular  element.  The 
validity  of  the  methods  developed  for  the  quadrilateral  shell  element  has  been  demonstrated 
by  numerical  results.  Specifically,  these  results  support  the  conclusions  drawn  from  previous 
work  that  rigid  body  modes  need  not  be  explicitly  included  in  the  element  displacement  functions 
if  the  latter  are  of  sufficiently  high  degree,  and  that  the  u,  v,  and  w  displacement  components 
should  be  of  equally  competent  degree.  Also  supported  by  the  results  is  the  validity  of  the 
method  adopted  for  approximating  the  shell  geometry. 

The  numerical  results  presented  include  solutions  to  both  static  and  dynamic  shell  re¬ 
sponse  problems.  These  have  led  to  the  following  conclusions  regarding  the  usage  of  the 
analysis  capability  developed. 

The  finite  element  static  results  show  excellent  agreement  with  theoretical  solutions. 
The  hemispherical  edge  effect  problem  illustrates  the  capability  of  the  doubly  curved  shell 
elements  to  determine  accurately  rapidly  varying  displacement  fields.  The  curved  elements 
also  predict  accurately  the  predominantly  membrane  behavior  for  the  internally  pressurized 
ellipse. 

The  dynamic  results  presented  show  the  effect  of  the  element-size-to-wave-length  ratio, 
and  of  the  different  ways  of  distributing  the  mass,  on  the  accuracy  of  the  analysis.  When  the 
generalized  mass  distribution  is  used,  the  size  of  the  eigenvalue  problem  may  be  reduced  by 
the  method  previously  described.  As  a  general  rule  it  appears  that  good  results  can  be  obtained 
by  reducing  up  to  90%  of  the  original  structural  freedoms,  retaining  only  the  normal  trans¬ 
lations,  w,  at  each  node.  Following  this  procedure,  with  an  element- size-to- wavelength  ratio 
of  1:4,  accuracy  of  better  than  5%  can  be  expected  in  the  predicted  frequency.  The  accuracy 
improves  rapidly  as  the  above  ratio  decreases,  but  may  also  deteriorate  rapidly  as  the 
ratio  increases.  The  latter  effect  was  noted  specifically  in  the  case  of  the  axisymmetric 
modes  of  the  60-degree  spherical  cap.  With  the  six-element  idealization  the  error  in  the  fre¬ 
quency  of  the  fourth  mode  obtained  with  no  reductions  was  9%;  however,  when  the  47  original 
degrees  of  freedom  were  reduced  to  the  11  normal  displacements  at  the  nodes  (an  80%  re¬ 
duction),  the  error  in  frequency  increased  to  115%.  Because  of  the  nature  of  the  mode  and 
the  placement  of  the  elements,  the  element- size-to- wavelength  ratio  is  greater  than  1:2  over 
a  large  region  near  the  edge  of  the  shell. 
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The  foregoing  statements  apply  to  cases  where  the  modal  response  is  of  the  typical  kind 
in  which  normal  displacements  participate  significantly.  In  cases  where  the  modes  of  interest 
are  of  a  type  in  which  normal  displacements  do  not  participate  significantly,  e.  g.,  in  torsional 
modes  of  a  cylinder,  then  u  or  v  displacement  components  should  be  retained  rather  than  w. 
Also,  it  is  presupposed  that  a  sufficient  number  of  elements  has  been  used  to  represent 
adequately  the  geometry  of  the  shell. 

The  lumped  mass  distribution,  while  not  elegant  in  the  mathematical  sense,  nevertheless 
has  certain  advantages: 

(1)  More  accurate  answers  are  obtained  with  a  coarser  idealization; 

(2)  Computing  time  necessary  to  generate  a  generalized  mass  matrix  is  saved: 

(3)  Reduction  of  the  mass  matrix,  which  is  roughly  four  times  slower  than  reduction  of  the 

stiffness  matrix,  is  unnecessary  and  considerable  computing  time  is  saved. 

For  this  reason  it  appears  that  the  use  of  lumped  masses  in  the  dynamic  analysis  of  shells 
is  more  attractive  from  the  engineering  viewpoint. 
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